#include "leafcurvature.h"

#include <iostream>
using namespace cv;

bool LeafCurvature::onBoundary(const Mat& leaf_mask,int x,int y)
{
    uchar min_value=leaf_mask.at<uchar>(x,y);
    uchar max_value=leaf_mask.at<uchar>(x,y);
    for(int i=-1;i<=1;i+=1){
        int xp=x+i;
        if(xp>=0 && xp<leaf_mask.rows){
            for(int j=-1;j<=1;j+=1){
                int yp=y+j;
                if(yp>=0 && yp<leaf_mask.cols){
                    uchar value=leaf_mask.at<uchar>(xp,yp);
                    min_value=min(min_value,value);
                    max_value=max(max_value,value);
                }
            }
        }
    }
    return min_value!=max_value;
}

ImageMask LeafCurvature::boundary(const Mat& leaf_mask)
{
    ImageMask boundary_mask(leaf_mask.size(),2);
    boundary_mask.mMask.setTo(0);

    for( int x = 0; x < leaf_mask.rows; x++ ) {
        for( int y = 0; y < leaf_mask.cols; y++ ) {
            if(leaf_mask.at<uchar>(x,y) && onBoundary(leaf_mask,x,y))
                boundary_mask.mMask.at<uchar> (x, y)=1;
        }
    }

    return boundary_mask;
}

ImageMask LeafCurvature::circleTemplate(float radius)
{
    int radius_int=ceil(radius);
    int template_size=radius_int*2+1;
    ImageMask circle_template(Size(template_size,template_size),2);

    Vec2f circle_center(radius_int,radius_int);
    for( int x = 0; x < circle_template.mMask.rows; x++ ) {
        for( int y = 0; y < circle_template.mMask.cols; y++ ) {
            bool some_inside=false;
            bool some_outsite=false;
            Vec2f cell_center(x,y);
            for(int i=-1;i<=1;i+=2){
                for(int j=-1;j<=1;j+=2){
                    Vec2f corner=cell_center+Vec2f(i,j)*0.5f;
                    if(norm(corner-circle_center)<radius)some_inside=true;
                    else some_outsite=true;
                }
            }
            if(some_inside&&some_outsite)
                circle_template.mMask.at<uchar>(x,y)=1;
            else
                circle_template.mMask.at<uchar>(x,y)=0;
        }
    }

    return circle_template;
}

ImageMask LeafCurvature::filledCircleTemplate(float radius)
{
    int radius_int=ceil(radius);
    int template_size=radius_int*2+1;
    ImageMask filled_circle_template(Size(template_size,template_size),2);

    Vec2f circle_center(radius_int,radius_int);
    for( int x = 0; x < filled_circle_template.mMask.rows; x++ ) {
        for( int y = 0; y < filled_circle_template.mMask.cols; y++ ) {
            Vec2f cell_center(x,y);
            if(norm(cell_center-circle_center)<radius)
                filled_circle_template.mMask.at<uchar>(x,y)=1;
            else
                filled_circle_template.mMask.at<uchar>(x,y)=0;
        }
    }

    return filled_circle_template;
}

std::vector<int> LeafCurvature::arclengthHistogram(const ImageMask& leaf_mask,float radius,int numberOfBins)
{
    std::vector<int> histogram;
    histogram.resize(numberOfBins);
    for(int i=0;i<numberOfBins;i++)histogram[i]=0;
//    Mat curvature(leaf_mask.size(),CV_8UC1);
//    curvature.setTo(0);

    ImageMask circle_template=circleTemplate(radius);
    int radius_int=(circle_template.mMask.rows-1)/2;

    int curvature_max=countNonZero(circle_template.mMask);
    int binStride=curvature_max/numberOfBins+1;

//    for( int x = 0; x < leaf_mask.rows; x++ ) {
//        for( int y = 0; y < leaf_mask.cols; y++ ) {
//            if(leaf_mask.at<uchar>(x,y) && onBoundary(leaf_mask,x,y))
//            {
    for(unsigned int p=0;p<leaf_mask.mBoundary.size();p++)
    {
        int x=leaf_mask.mBoundary[p].y;
        int y=leaf_mask.mBoundary[p].x;
                int count=0;
                int start_x=x-radius_int;
                int start_y=y-radius_int;
                for( int i = 0; i < circle_template.mMask.rows; i++ ) {
                    for( int j = 0; j < circle_template.mMask.cols; j++ ) {
                        int image_i=start_x+i;
                        int image_j=start_y+j;
                        if(circle_template.mMask.at<uchar>(i,j)>0 &&
                                image_i>=0 && image_i<leaf_mask.mMask.rows &&
                                image_j>=0 && image_j<leaf_mask.mMask.cols &&
                                leaf_mask.mMask.at<uchar>(image_i,image_j)>0){
                            count++;
                        }
                    }
                }
                histogram[count/binStride]++;
                //curvature.at<uchar>(x,y)=count;
    }
//            }
//        }
//    }

    //std::cout<<curvature_max<<std::endl;
    //double min_value,max_value;
    //minMaxLoc(curvature,&min_value,&max_value);
    //vstd::cout<<min_value<<" "<<max_value<<std::endl;

//    for( int x = 0; x < curvature.rows; x++ ) {
//        for( int y = 0; y < curvature.cols; y++ ) {
//            if(leaf_mask.at<uchar>(x,y) && onBoundary(leaf_mask,x,y))
//                histogram[curvature.at<uchar>(x,y)/binStride]++;
//        }
//    }

    return histogram;
}

std::vector<int> LeafCurvature::areaHistogram(const ImageMask& leaf_mask,float radius,int numberOfBins)
{
    std::vector<int> histogram;
    histogram.resize(numberOfBins);
    for(int i=0;i<numberOfBins;i++)histogram[i]=0;
//    Mat curvature(leaf_mask.mMask.size(),CV_8UC1);
//    curvature.setTo(0);

    ImageMask filled_circle_template=filledCircleTemplate(radius);
    int radius_int=(filled_circle_template.mMask.rows-1)/2;

    int curvature_max=countNonZero(filled_circle_template.mMask);
    int binStride=curvature_max/numberOfBins+1;

//    for( int x = 0; x < leaf_mask.rows; x++ ) {
//        for( int y = 0; y < leaf_mask.cols; y++ ) {
//            if(leaf_mask.at<uchar>(x,y) && onBoundary(leaf_mask,x,y))
//            {
    for(unsigned int p=0;p<leaf_mask.mBoundary.size();p++)
    {
        int x=leaf_mask.mBoundary[p].y;
        int y=leaf_mask.mBoundary[p].x;
                int count=0;
                int start_x=x-radius_int;
                int start_y=y-radius_int;
                for( int i = 0; i < filled_circle_template.mMask.rows; i++ ) {
                    for( int j = 0; j < filled_circle_template.mMask.cols; j++ ) {
                        int image_i=start_x+i;
                        int image_j=start_y+j;
                        if(filled_circle_template.mMask.at<uchar>(i,j)>0 &&
                                image_i>=0 && image_i<leaf_mask.mMask.rows &&
                                image_j>=0 && image_j<leaf_mask.mMask.cols &&
                                leaf_mask.mMask.at<uchar>(image_i,image_j)>0){
                            count++;
                        }
                    }
                }

                histogram[count/binStride]++;
                //curvature.at<uchar>(x,y)=count;
    }
//            }
//        }
//    }

//    for( int x = 0; x < curvature.rows; x++ ) {
//        for( int y = 0; y < curvature.cols; y++ ) {
//            if(leaf_mask.at<uchar>(x,y) && onBoundary(leaf_mask,x,y))
//                histogram[curvature.at<uchar>(x,y)/binStride]++;
//        }
//    }

    return histogram;
}

std::vector<float> LeafCurvature::normalizedArclengthHistogram(const ImageMask& leaf_mask,float radius,int numberOfBins)
{
    std::vector<int> original=arclengthHistogram(leaf_mask, radius, numberOfBins);
    int sum=0;
    for(unsigned int i=0;i<original.size();i++)sum+=original[i];
    std::vector<float> output(original.size());
    for(unsigned int i=0;i<original.size();i++)output[i]=(float)original[i]/(float)sum;
    return output;
}

std::vector<float> LeafCurvature::normalizedAreaHistogram(const ImageMask& leaf_mask,float radius,int numberOfBins)
{
    std::vector<int> original=areaHistogram(leaf_mask, radius, numberOfBins);
    int sum=0;
    for(unsigned int i=0;i<original.size();i++)sum+=original[i];
    std::vector<float> output(original.size());
    for(unsigned int i=0;i<original.size();i++)output[i]=(float)original[i]/(float)sum;
    return output;
}
